Información georeferenciada en repositorios online

Vamos a buscar data de repositorios online, sobretodo usaremos Open Street Map (OSM)

# Cargamos las libraries

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(ggmap)
## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
##   Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
##   OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

OSM tiene informacion de todo el mundo y vamos a querer bajar solo un área de interés. Vamos a ponerle como string la query de la búsqueda del lugar.

La wiki de OSM tiene la info necesaria para entender el catálogo: https://wiki.openstreetmap.org/wiki/ES:Objetos_del_mapa

# Obtenemos la bounding box del área de interes

bbox_aoi <- getbb("Comuna 10, Ciudad Autonoma de Buenos Aires, Buenos Aires, Argentina")
bbox_aoi
##         min       max
## x -58.53145 -58.47131
## y -34.64581 -34.60912

Chequeemos que las coordenadas de la bounding box estan ok haciendo un mapa y viendo que sea la zona del mundo que nos interesa, no sea cosa que estemos en otro hemisferio porejemplo. Podría pasar, se imaginan?

register_stadiamaps("cebdf63b-1fbd-40ac-af86-fc993f352022", write = TRUE)
## ℹ Replacing old key (cebdf63b) with new key in /Users/danielarisaro/.Renviron
mapa_aoi <- get_stadiamap(bbox = bbox_aoi,
                              maptype = "stamen_toner_lite",
                              zoom = 14)
## ℹ © Stadia Maps © Stamen Design © OpenMapTiles © OpenStreetMap contributors.
ggmap(mapa_aoi)

Ademas del area de la comuna tambien podemos usar los poligonos que las delimitan

poligono_aoi <- getbb("Comuna 10, Ciudad Autonoma de Buenos Aires, Buenos Aires, Argentina",
                          format_out = "sf_polygon")

Podemos hacer la superposicion del area y del poligono para ver como lucen juntos

ggmap(mapa_aoi) +
  geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 1.25, inherit.aes = FALSE) +
  labs(title = "Comuna 10",
       subtitle = "Ciudad Autonoma de Buenos Aires",
       caption = "Fuente: OpenStreetMap.") +
  theme_void()
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

How to get information from OSM?

La informacion disponible de OSM la podemos consultar con dos funciones: available_features() y available_tags() Las features nos dan keys y cada key tiene tags

all_features <- available_features()

print(all_features[1:10])
##  [1] "4wd_only"                "abandoned"              
##  [3] "abutters"                "access"                 
##  [5] "addr"                    "addr:city"              
##  [7] "addr:conscriptionnumber" "addr:country"           
##  [9] "addr:county"             "addr:district"
available_tags('building')
## # A tibble: 101 × 2
##    Key      Value          
##    <chr>    <chr>          
##  1 building allotment_house
##  2 building annexe         
##  3 building apartments     
##  4 building bakehouse      
##  5 building barn           
##  6 building barracks       
##  7 building beach_hut      
##  8 building boathouse      
##  9 building bridge         
## 10 building bungalow       
## # ℹ 91 more rows

Como opciones podemos descargar lineas, puntos o poligonos (y multipolígonos). Vamos de a uno?

Puntos!

Vamos a descargar algunos puntos en particular de las churchs

churchs_aoi <- opq(bbox_aoi)
churchs_aoi <- add_osm_feature(churchs_aoi,
                                       key = "building",
                                       value = c("building", "church"))

churchs_aoi <- osmdata_sf(churchs_aoi)
churchs_aoi
## Object of class 'osmdata' with:
##                  $bbox : -34.6458134,-58.5314494,-34.609115,-58.4713084
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 126 points
##             $osm_lines : NULL
##          $osm_polygons : 'sf' Simple Features Collection with 12 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL

Me quedo con los puntos. Me podria quedar con los poligonos tambien haciendo osm_polygons

churchs_aoi <- churchs_aoi$osm_points

A ver cuantos puntos tenemos:

dim(churchs_aoi)
## [1] 126   3

126 puntos!

ggmap(mapa_aoi) +
  geom_sf(data = churchs_aoi, inherit.aes = FALSE) +
  geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
  labs(title = "Edificios religiosos - iglesias",
       subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
       caption = "Fuente: OpenStreetMap") +
  theme_void() +
  theme(title = element_text(size = 8, face = "bold"), #tamaño de titulo del mapa
        plot.caption = element_text(face = "italic", colour = "gray35", size = 7)) #tamaño de nota al pie
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Vamos a ver exactamente cuantas churchs entran en la comuna 10. Para eso vamos a usar la st_intersection()

churchs_aoi <- st_intersection(churchs_aoi, poligono_aoi)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
dim(churchs_aoi)
## [1] 68  3

son sesentayocho. che me parece UN MONTON. Y por qué hay puntos tan contiguos en el mapa? Pues no lo se, pero habría que ver qué entiende OSM por church

Lineas, lineas lineas,

Tambien podemos descargar lineas del area of interest.

calles_aoi <- opq(bbox_aoi)

calles_aoi <- add_osm_feature(calles_aoi,
                                  key = "highway",
                                  value = c("motorway", "trunk", "primary", "secondary", "tertiary"))

calles_aoi <- osmdata_sf(calles_aoi)
calles_aoi
## Object of class 'osmdata' with:
##                  $bbox : -34.6458134,-58.5314494,-34.609115,-58.4713084
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 2701 points
##             $osm_lines : 'sf' Simple Features Collection with 718 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 0 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL

Nos quedamos con las lines

calles_aoi <- calles_aoi$osm_lines
ggmap(mapa_aoi) +
  geom_sf(data = calles_aoi, color = "#2b2d42", inherit.aes = FALSE, lwd = 0.7) +
  geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
  labs(title = "Calles principales",
       subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
       caption = "Fuente: OpenStreetMap") +
  theme_void() +
  theme(title = element_text(size = 8, face = "bold"), #tamaño de titulo del mapa
        plot.caption = element_text(face = "italic", colour = "gray35", size = 7)) #tamaño de nota al pie
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Por último, polígonos.

Para finalizar vamos a hacer una query de poligonos (overpass query, remember). Acá me anoto una tarea que es ver bien el tema de los tags porque no los entiendo del todo: en los values estan concatenadas diferentes keys o hay una adentro de la otra?

parques_aoi <- opq(bbox_aoi)

parques_aoi <- add_osm_feature(parques_aoi,
                                   key = "leisure",
                                   value = c("park", "garden", "nature_reserve"))

parques_aoi <- osmdata_sf(parques_aoi)
parques_aoi
## Object of class 'osmdata' with:
##                  $bbox : -34.6458134,-58.5314494,-34.609115,-58.4713084
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 1879 points
##             $osm_lines : 'sf' Simple Features Collection with 48 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 110 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : 'sf' Simple Features Collection with 8 multipolygons

Ahora vamos a tomar los poligonos y multipoligonos de este objeto

parques_polygons <- parques_aoi$osm_polygons
parques_multipolygons <- parques_aoi$osm_multipolygons

dim(parques_polygons)
## [1] 110  29
ggmap(mapa_aoi) +
  geom_sf(data = parques_polygons, fill = "#588157", color = NA, inherit.aes = FALSE) +
  geom_sf(data = parques_multipolygons, fill = "#588157", color = NA, inherit.aes = FALSE) +
  geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
  labs(title = "Espacios verdes",
       subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
       caption = "Fuente: OpenStreetMap") +
  theme_void()
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Tenemos dos elementos que son los poligonos y los multipoligonos y queremos que sea una sola! Para eso vamos a seleccionar by the id and then we will bind them together oh yeah

parques_polygons <- parques_polygons %>%
  select(osm_id, name)
parques_multipolygons <- parques_multipolygons %>%
  select(osm_id, name)
parques_aoi <- rbind(parques_polygons, parques_multipolygons)

Ahora intersectamos el poligono de la comuna de interes con los poligonos de los parques verdes

parques_aoi <- st_intersection(parques_aoi, poligono_aoi)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
ggmap(mapa_aoi) +
  geom_sf(data = parques_aoi, fill = "#588157", color  =NA, inherit.aes = FALSE) +
  geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
  labs(title = "Espacios verdes",
       subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
       caption = "Fuente: OpenStreetMap") +
  theme_void() +
  theme(title = element_text(size = 8, face = "bold"), #tamaño de titulo del mapa
        plot.caption = element_text(face = "italic", colour = "gray35", size = 7)) #tamaño de nota al pie
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

jajajaja me quedan los parques on top de los names pero ya fue

Vamos a calcular la superficie de las áreas verdes

parques_aoi <- parques_aoi %>%
  mutate(area_m2 = st_area(parques_aoi))

summary(parques_aoi$area_m2)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    10.55   239.33   876.10  2604.68  2479.92 18665.13

El area promedio de los espacios verdes en la comuna 10 es de 2600 m2 creo yo que es muy poco poquisimo, terrible mala performance

summarise(parques_aoi, area_m2 = sum(area_m2))
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.52973 ymin: -34.64361 xmax: -58.47425 ymax: -34.614
## Geodetic CRS:  WGS 84
##          area_m2                       geometry
## 1 210979.2 [m^2] MULTIPOLYGON (((-58.4916 -3...

El área total de los parques es de 210979 m2, 0.2 km2